The phase diagram of ice Ih, II, and III: a quasi-harmonic study 
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The phase diagram of ice Ih, II, and III is studied by a quasi-harmonic approximation. The results 
of this approach are compared to phase diagrams previously derived by thermodynamic integration 
using path integral and classical simulations, as well as to experimental data. The studied models 
are based on both flexible (q-TIP4P/F) and rigid (TIP4P/2005, TIP4PQ/2005) descriptions of the 
water molecule. Many aspects of the simulated phase diagrams are reasonably reproduced by the 
quasi-harmonic approximation. Advantages of this simple approach are that it is free from the 
statistical errors inherent to computer simulations, both classical and quantum limits are easily 
accessible, and the error of the approximation is expected to decrease in the zero temperature 
limit. We find that the calculated phase diagram of ice Ih, II, and III depends strongly on the 
hydrogen disorder of ice III, at least for cell sizes typically used in phase coexistence simulations. 
Either ice II (in the classical limit) or ice III (in the quantum one) may become unstable depending 
upon the proton disorder in ice III. The comparison of quantum and classical limits shows that 
the stabilization of ice II is the most important quantum effect in the phase diagram. The lower 
vibrational zero-point energy of ice II, compared to either ice Ih or III, is the microscopic origin 
of this stabilization. The necessity of performing an average of the lattice energy over the proton 
disorder of ice III is discussed. 

PACS numbers: 64.60.-i,64.60.De, 63.20.-e, 63.20.Ry 



I. INTRODUCTION 

An outstanding property of water is the diversity of 
ice phases that are found in its phase diagram^ Six- 
teen different crystalline ice phases have been identified 
so far, a number that is likely to increase in the fu- 
ture. In all phases, except ice X, the water molecule 
appears as a well defined entity that is part of a network 
of molecules connected by H-bonds. In this network each 
water molecule is surrounded by four others in a more or 
less distorted tetrahedral coordination. The orientation 
of each molecule with respect to its four nearest neighbors 
satisfies the Bernal-Fowler ice rules. They state that the 
H2O molecule is oriented so that its two protons point 
toward adjacent oxygen atoms and that there must be 
one and only one proton between two adjacent oxygen 
atoms' 

The existence of orientational disorder in the water 
molecules is a property of several ice phases. While the 
oxygen atoms display a full occupancy (/) of their crys- 
tallographic positions, the hydrogen atoms may display 
a disordered spatial distribution as evidenced by a frac- 
tional occupancy of their lattice sites. Ice Ih, the stable 
phase of ice under normal conditions, displays full proton 
disorder compatible with the Bernal-Fowler rules, i.e., 
occupancies of H-sites of / = 0.5. However ice II is H- 
ordered, while ice III is characterized by a partial proton 
ordering, i.e., some fractional occupancies of H-sites are 
different from 0.5. Order-disorder transitions have been 
observed for several pairs of ice phases (Ih-XI, III-IX, 
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V-XIII, VI-XV, VII- VIII, XII-XIV). The orientational 
ordering implies a whole reorganization of the H-bond 
network that is kinetically unfavorable. In several disor- 
dered phases this transition only occurs after doping with 
either bases (in the case of ice Ih) or acids (for ices V, VI, 
and XII) . The creation of defects provides a mechanism 
favoring the rearrangement of the H-bond network. 3 

The simulation of the complex phase diagram of 
water is an interesting challenge. Large portions of 
the phase diagram have been simulated using rigid 
models (TIP4P, 4 TIP4P/2005, 5 and TIP4PQ/2005 6 ), 
and smaller regions using a flexible water model (q- 
TIP4P/F)!. Let us present a brief summary of these 
TIP4P-like models. The TIP4P potential is based on a 
point charge description of a rigid water molecule supple- 
mented by an additional Lennard- Jones interaction be- 
tween the oxygen centers. It was parameterized by Jor- 
gensen et al. in 1983. 8 An optimized parameterization 
of the same model was labeled as TIP4P/2005. 5 Both 
model potentials have been employed in classical simula- 
tions. Consideration of quantum effects by path integral 
simulations with the TIP4P/2005 model led to unphys- 
ical results, e.g., ice II was predicted to be more stable 
than ice Ih at low temperatures. Then, a small increase in 
the point charges was proposed to avoid this problem and 
the new parameterization was labeled as TIP4PQ/2005. 9 
An interesting recent contribution was to add to the rigid 
TIP4P/2005 model an anharmonic potential energy term 
to treat the molecular flexibility of water in quantum 
simulations, giving rise to the q-TIP4P/F models Sum- 
marizing, TIP4P/2005 and TIP4PQ/2005 are rigid wa- 
ter models intended for its use in classical and quantum 
simulations, respectively, while q-TIP4P/F is a flexible 
model for quantum simulations. 
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A comprehensive review of the calculation of free ener- 
gies in water phases with the thermodynamic integration 
(TI) method can be found in Ref. [HI The classical phase 
diagram of water, simulated with the TIP4P/2005 model, 
shows a reasonable qualitative agreement to the experi- 
mental one, in particular in the complex region of stabil- 
ity of ices Ih, II, III, V, and VI. - The phase diagram of ice 
Ih, II, and III has been also investigated using the flexible 
water model (q-TIP4P/F) in the classical limits It was 
found that ice II is unstable, as its stability region was oc- 
cupied by ice III. A plausible explanation of the difference 
in the stability of ice II found with the TIP4P/2005 and 
q-TIP4P/F models might be that the geometry, dipole 
and quadrupole moments of the molecules in the flexible 
model can vary between the different ice phases, which 
is not the case for a rigid model. The simulation of the 
phase diagram of water using quantum path integral sim- 
ulations with the TIP4PQ/2005 model has been recently 
reported for a temperature range between 140 and 300 K, 
and pressures up to 1.2 GPa£ The quantum results were 
compared to the classical expectation. One of the most 
streaking difference between the classical and quantum 
limits is the region of stability of ice III, which is much 
lower (and in better agreement to the experiment) in the 
quantum case. 

Despite the overall agreement found between simulated 
and experimental phase diagrams of water, it is obvious 
that some properties can not be reproduced by the em- 
ployed empirical models. Order-disorder transitions are 
a prominent example. It is well documented that these 
water models are unable to predict that ice XI (the ferro- 
electric ordered counterpart of ice Ih with Cmc2i spatial 
symmetry) is the most stable ice phase below 72 K at at- 
mospheric pressured Also the order-disorder transition 
between ice VII and VIII is poorly reproduced by the 
empirical models^ It seems that these effective potentials 
fail to describe the energetics of proton rearrangements in 
ice. Therefore the location of order-disorder transitions 
and the identity of the ordered low temperature phases 
is inadequately predicted. 3 An interesting question is to 
what extent the use of ab initio density functional the- 
ory (DFT) can improve these limitations. In this respect, 
DFT studies of liquid water and ice have revealed seri- 
ous differences with experimental data in both diffusive 
and structural properties that seems to be related to the 
subtle contribution of van der Waals dispersion forces to 
the cohesive energy of the water phases. The application 
of new functionals specially designed to treat van der 
Waals interactions is focus of recent interest in modeling 
water 3^ 

The quasi-harmonic (QH) approximation (QHA) al- 
lows to compute the partition function of a solid phase 
as an analytic function of the crystal volume and the 
temperature; 1 ® Some advantages of this approach are the 
straightforward derivation of equilibrium thermodynamic 
properties, the absence of statistical errors (as opposed 
to any simulation method) and the possibility to account 
for finite size effects by a Brillouin zone integration of the 



phonon dispersion curves, rather than by increasing the 
size of the cell. The QHA in combination with ab ini- 
tio DFT has allowed the explanation of the inverse iso- 
tope effect in the crystal volume of ice Ih at atmospheric 
pressured Also the negative thermal expansion of ice 
Ih at low temperatures has been studied by the QHA, 18 
as well as the elastic moduli and mechanical stability of 
the H-ordered ice VIII/^ In addition, the mechanical sta- 
bility of ice Ih under pressure has been studied by this 
approximation^ The ice VII- VIII phase boundary has 
been studied by a QHA in a 16-molecule supercell with 
ab initio DFT calculations of total energies and phonon 
frequencies i^i 

The validity of the QHA is restricted by the possible 
presence of anharmonic effects beyond those included in 
the approximation. Thus, a direct check of the QHA is 
the comparison to numerical simulations that fully con- 
sider the anharmonicity of the interatomic interactions. 
The QHA prediction of the volume, enthalpy, kinetic en- 
ergy, and heat capacity, of ice Ih, II, and III has been 
compared to both quantum and classical simulations us- 
ing the q-TIP4P /F models The comparison in a (T, P) 
range up to 300 K and 1 GPa shows a remarkable over- 
all agreement for the three ice phases. An interesting 
aspect of the QHA is that it is sensible enough to pre- 
dict differences in the anharmonicity of different water 
models, that are in agreement to the corresponding fully 
anharmonic simulations. For example, the QHA predicts 
that the thermal expansion of ice Ih at low temperatures 
is negative for the q-TIP4P/F and TIP4P models but 
positive (or slightly negative) for the TIP5P and ST2 
potentials; 22 ! 23 Moreover, the isotope effect in the crystal 
volume of ice Ih is predicted by the QHA to be anomalous 
(as in the experiment) with a DFT functional, but nor- 
mal with the q-TIP4P/F model^ We stress that these 
QHA predictions of anharmonic effects are in agreement 
to results of available computer simulations. It may be 
somewhat surprising that the simple QHA approxima- 
tion is able to reproduce the anharmonicity of complex 
ice phases with a similar accuracy as that shown for solids 
with much simpler crystal structures such as noble gases 
and elemental semiconductors (Si, Ge)^ — 

The purpose of the present paper is to check the ca- 
pability of the QHA to predict the phase diagram of ice 
Ih, II, and III. The layout of the manuscript is as fol- 
lows. A summary of the employed computational con- 
ditions is presented in Sec. [TTJ The generation of the 
ice structures is introduced in Sec. IIIII The QH phase 
diagram is studied for the flexible q-TIP4P/F model in 
Sec. HVl The results for the rigid models, TIP4P/2005 
and TIP4PQ/2005, are presented in Sees. [V] and ED re- 
spectively. Our main focus of interest is the influence 
of proton disorder in the calculated phase diagram, the 
comparison of the QHA to previous simulation results, 
and the differences between the quantum and classical 
limits. The necessity of performing disorder averaging is 
discussed in Sec. I VIII The paper closes with the conclu- 
sions. 
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II. COMPUTATIONAL CONDITIONS 

The QHA employed for the ice phases has been in- 
troduced in Ref. [22- We present here a brief summary. 
The Helmholtz free energy of an ice phase with N water 
molecules in a cell of volume V and at temperature T is 
given by 

F(V,T) = U S (V) + F V (V,T)-TS H , (1) 

where Us(V) is the static zero-temperature classical en- 
ergy, i.e., the minimum of the potential energy when the 
volume of the ice cell is V. The entropy Sh is related to 
the disorder of hydrogen and it vanishes for ordered ice 
phases as ice II. Sh was estimated by Pauling for fully 
disordered phases as 2 ^ 

S H =Nk B \n^. (2) 

F V (V,T) is the vibrational contribution to F. In the 
quantum limit is given by 

F t ,(t/,T) = ^^ + iln[l-exp(-^)]) ■ (3) 

Here j3 is the inverse temperature: l/fc^T. uj^ are the 
wavenumbers of the harmonic lattice vibrations for the 
volume V, with h combining the phonon branch index 
and the wave vector within the Brillouin zone. In the 
classical limit the vibrational contribution amounts to 

F VtCla (V,T)=J2l^(^k) ■ (4) 

k P 

The Gibbs free energy, G(T, P), is obtained by seeking for 
the volume, V m i ni that minimizes the function F(V, T) + 
PV, as 

G(T, P) = F(V mm ,T) + PV mm . (5) 

The implementation of the QHA for an ice phase follows 
these steps j 2 ^ 

i) Find the reference cell that minimizes the static en- 
ergy Us- This minimization implies optimization of both 
cell shape and atomic positions. The resulting volume is 
V re f and the corresponding static energy Us,ref- 

ii) Select a grid of 50 volumes in a range of inter- 
est [V m in,V max }. The ice cell with volume Vi is set by 
isotropic scaling of the reference cell. Subsequently, each 
ice cell is held fixed while minimizing the static energy 
Us(Vi) with respect to the atomic positions. The crystal 
phonons, u>k(Vi), are obtained after the minimization. 

Hi) Calculate the function F(Vi,T) by Eq. (JTJ. The 
minimum of F(Vi, T) as a function of V is determined by 
a fit to a 5th degree polynomial in V. 

The phase diagram of the ice phases is derived by a 
brute force method, i.e., given a state point (T,P) one 
calculates the Gibbs free energy of all ice phases and then 



the stable phase is selected as the one with the lowest 
value of G. 

The crystal phonon calculation has been performed 
by the small-displacement methodi 28 ' 29 For the flexible 
water model the atomic displacement employed in this 
work is Sx = 10~ 6 A along each Cartesian direction. 
For the rigid water models the molecular displacements 
imply translations (by 10~ 6 A along the Cartesian di- 
rections) and rotations (by 10~ 7 rad along the Carte- 
sian axes) of the rigid molecules. See Ref. [30] for a 
full account of the calculation of the external phonon 
modes associated to rigid units. We have used a T sam- 
pling (k = 0) of the crystal phonons, as this condition is 
implicitly assumed in simulation studies using periodic 
boundary conditions^ The Lennard- Jones interaction 
between oxygen centers was truncated at a distance of 
r c = 8.5 A, and standard long-range corrections for both 
potential energy and pressure were computed assuming 
that the pair-correlation function is unity for r > r c . 31 
Long-range electrostatic potential and forces were calcu- 
lated with the Ewald method. 

As a check for the assumption of isotropic expansion of 
the reference cell made in the step »), we have performed 
a QHA of the Gibbs free energy of ice II by relaxing 
this constraint. To this aim we have derived a set of 50 
cell volumes V, by performing an optimization of both 
cell shapes and atomic positions at 50 different pressures 
in the range [-1.7, 3] GPa. The new result for the free 
energy of ice II reveals that the assumption of isotropic 
scaling of the reference cell modifies the values of G by 
less than 0.01 kJ/mol. This change in G has a small 
effect in the phase diagram of ice reflected by rigid shifts 
of the calculated coexistence lines of ice II by less than 
2 K. 



III. ICE STRUCTURES 

Supercells of similar size to those employed in recent 
simulations 6,7 have been used in the QH derivation of 
the phase diagram. The number of molecules were N= 
288 for ice Ih, and N= 324 for ice II and III. The ice Ih 
cell was orthorhombic with parameters (4ai , 3\/3ai , 3a3), 
with (ai,a3) being the standard hexagonal lattice vec- 
tors of ice Ihj22 Ice II and III were studied by 3 x 3 x 3 
supercells of the crystallographic cell, which belong to 
the rhombohedral and the tetragonal crystal systems, 
respectively] 33 i 34 While ice II is proton ordered both ice 
Ih and III display orientational disorder of the water 
molecules. The algorithm proposed by Buch et al. was 
applied for the random generation of full proton disor- 
dered structures (/ = 0.5) with vanishing cell dipole 
moment fi 2 " In the case of ice III, the neutron diffrac- 
tion experiments show that only 1 /3 of the H-sites has 
/ = 0.5, while the other 2/3 display occupancies of 
/ = 0.35 and / = 0.65, respectively^ The Buch's al- 
gorithm has been modified for the generation of random 
structures with partial H-disorder, i.e., having fractional 
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Figure 1: Minimum static energy and corresponding cell vol- 
ume of randomly generated H-isomers of ice III. The H- 
isomers display either full (open circles) or partial (closed 
squares) H-disorder. The results correspond to the q- 
TIP4P/F model for a 324-molecule supercell. An arrow points 
to the H-isomer with full H-disorder and largest energy. This 
particular isomer is employed in the calculation of phase dia- 
grams with different water models. The line is a linear fit to 
the data. 



occupancies different from / = 0.5^ 

An interesting practical question is the importance of 
proton disorder in the evaluation of the partition function 
of the H-disordered phases, using a single H-isomer. This 
point has not been addressed earlier in the simulation of 
phase diagrams of ice using either TIP4P, TIP4P/2005, 
TIP4PQ/2005, or q-TIP4P/F models^! These simula- 
tions used a single H-isomer of the disordered phases (ice 
Ih and III), generated by a random algorithm. Moreover, 
the H-isomer for ice III was selected either with partial 
H-disorder in the earlier simulation with the rigid TIP4P 
model 4 or with full H-disorder in more recent simulations 
with the q-TIP4P/F and TIP4PQ/2005 potentials^ 
The question to be addressed here is how large might be 
the influence of the selected H-isomer in the calculated 
phase diagram. 

To this aim we have generated a random set of six 
H-isomers with full H-disorder and vanishing cell dipole 
moment for ice III. The result of their energy minimiza- 
tion with the flexible q-TIP4P /F potential is represented 
in Fig. [T] The volume, V re f, and the corresponding 
minimized potential energy, Us, r ef; is displayed by open 
circles for each H-isomer. We note that the static energy, 
Us,ref, of the six H-isomers spreads in an energy window 
of about 0.3 kJ/mol. The volume, V^ e /, and the min- 
imized potential energy, Us,ref, are related in a nearly 
linear way. A second set of six random isomers with 
vanishing cell dipole moment has been generated by im- 
posing the partial H-disorder encountered in the diffrac- 
tion experiment of ice III. 34 The results of the energy 



minimization for the partially disordered structures are 
presented as closed squares in Fig. [TJ All isomers having 
partial H-disorder display larger static energy than the 
isomers with full H-disorder. 

Similar behavior to that shown in Fig. [TJis found if the 
minimization of the energy of the H-isomers is performed 
with the rigid TIP4P/2005 model. The main difference is 
that the dispersion of the Us,ref values increases slightly. 
Thus, two main conclusions can be derived from the re- 
sults of the energy minimization in Fig. [TJ The first is 
that the energetics associated to full versus partial H- 
disorder in ice III is incorrectly described by effective 
TIP4P-likc models, i.e., full disorder is predicted to be 
more stable than partial one. Note that the configura- 
tional entropy Sr will help to stabilize further the full 
disordered ice at any finite temperature, since Sh in this 
case is larger than for the partial H-disorder. This be- 
havior is in contradiction to the partial H-disorder exper- 
imentally found for ice Illi^i Our result is in line with the 
reported limitations of these effective potentials to repro- 
duce the energetics of the H-bond rearrangement in the 
order-disorder transition of ice Ih-XI and VH-VIIIi&i^ 

Our second conclusion is that the large dispersion in 
Us,ref obtained for ice III using cells with 324 molecules 
must affect the phase diagram whenever it is calculated 
with a single random H-isomer. The energy dispersion 
is caused by the proton disorder in the H-isomers. The 
sampling of this (large) dispersion of static cell energies 
with a single H-isomer can be considered the origin of a 
finite size effect. In the thermodynamic limit, the energy 
distribution of Us, re f should approximate a delta func- 
tion centered at the energy corresponding to the most 
probable H-bond distribution. Therefore, the finite size 
effect caused by the insufficient sampling of the proton 
disorder (or the Us.ref energies) with a single H-isomer is 
expected to decrease as the size of the cell increases. An 
alternative way to reduce this finite size effect is to make 
a disorder averaging of the lattice energy of the employed 
ice cell. This point will be commented in Sec. IVIII 

In the case of ice Ih we have also generated a ran- 
dom set of six H-disordered structures. However, in this 
case the static energy, Us,refi of the H-isomers varies 
in a rather small energy window (of about 0.01 kJ/mol 
for the q-TIP4P/F model) and the corresponding volume 
changes by less than 0.04%. Thus, finite size effects re- 
lated to the H-disorder are expected to be much lower in 
ice Ih than in ice III. 

A comparison of V re f and Us,ref calculated with the 
q-TIP4P/F model for ice Ih, II, and III are given in Tab. 
UJ The data for ice III correspond to the H-isomer labeled 
with an arrow in Fig. [TJ The classical internal energy at 
zero temperature and pressure (T = 0, P = 0) is Us, r ef- 
The QH result in the quantum limit for the ice volume 
(Vb), static energy {Us,o), and zero-point energy (Uz.o) 
at T — and P = are also summarized in Tab. [TJ The 
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Table I: Volume (V re f) and static energy (Us,ref) of the minimum energy configuration of the studied ice phases. The quantum 
QH results for the volume (Vb), static energy (Us,o), zero-point energy (Uz,o), and internal energy (Uo) are also given at T = 
and P — 0. The data for ice III correspond to the H-isomer labeled by an arrow in Fig. [1] All results were derived with the 
q-TIP4P/F model. The last two columns show the difference with the data of ice II. [Vmin, V m ax] is the volume interval studied 
by the QHA for each phase. 



X (q-TIP4P/F) 


Ih 


II 


III 


AX (Ih-II) AX (III-II) 


V re f (A 3 /molec.) 


30.96 


24.14 


24.99 


6.82 


0.85 


U s , re f (kJ/mol) 


-61.98 


-60.84 


-60.86 


-1.14 


-0.02 


V (A 3 /molec.) 


32.23 


25.11 


25.90 


7.12 


0.79 


Us, (kJ/mol) 


-61.74 


-60.60 


-60.77 


-1.14 


-0.17 


Uz,o (kJ/mol) 


68.75 


68.08 


68.72 


0.67 


0.64 


Uo (kJ/mol) 


7.01 


7.47 


7.95 


-0.46 


0.47 


Vmin (A 3 /molec.) 
Vmax (A 3 /molec.) 


29.47 


21.75 


22.48 






35.05 


27.31 


28.22 







zero-point energy, Uz.o, is calculated as 

hu k (V ) 



(G) 



Note that the zero-point energy of ice II is lower than 
that of ice Ih and III. In the quantum limit, the internal 
energy of the ice phases at T — and P = is 



U = U S fl + Uz,o 



(7) 



The ice structures studied in Tab. [I] have been analyzed 
also with the rigid TIP4P/2005 and TIP4PQ/2005 mod- 
els. The corresponding results are presented in Tabs. |n] 
and Hill Note that the zero-point energy (Uz,o) of the 
rigid models is about four times smaller than that of 
flexible water because of the neglect of intramolecular 
motion. 



IV. FLEXIBLE Q-TIP4P/F MODEL 

In this section the QH phase diagram of ice Ih, II, 
and III is derived with the q-TIP4P/F model in both 
classical and quantum limits. Studied temperatures are 
in the interval [0,300 K] and pressures in the range [0, 
0.35 GPa]. The calculation is done for each of the six 
random H-isomers of ice III having full H-disorder. A 
comparison to available TI simulations is provided in the 
classical limit. 7 



A. Classical limit 

The QH phase diagram of ice Ih, II, and III calculated 
in the classical limit is plotted in Fig. [5J Coexistence 
lines are displayed for each of the six studied H-isomers of 
ice III as continuous curves. We find that the finite size 
effect related to the H-disorder in ice III is important. 
In particular, the area where ice III is stable strongly 



0.1 



~l — i — i — i — i — | — i — i — i — r~ 
III (6 H-isomers) 




cla. QHA 

cla. TI (Habershon et al.) 
TP (Ih-II-ffl) 



q-TIP4P/F 



100 200 
T(K) 



300 



Figure 2: Phase diagram of ice Ih-II-III calculated with the 
q-TIP4P/F model in the classical limit. The full lines show 
the QH results. The QH calculation is done for six randomly 
chosen H-isomers having full H-disorder in a cell with 324 
water molecules. The dotted lines are the classical TI results 
of Ref. that include also the boundary with the liquid (L) 
phase. Circles show the position of the triple point (TP) for 
ice Ih-II-III. 



depends upon the H-isomer. The dotted lines in Fig. [5] 
show the coexistence lines for ice Ih-III, Ih-liquid, and 
Ill-liquid as derived from the classical TI simulations of 
Ref. @with the q-TIP4P/F model. The coexistence line 
Ih-III is parallel to our QHA results. 

For the various H-isomers, the ice Ih-III phase bound- 
ary is shifted by a nearly constant pressure. The disper- 
sion of the static energy, Us. ref, is the factor responsible 
for the different phase behavior of the H-isomers. Vibra- 
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Table II: Volume and energies of the studied ice phases as derived with the rigid TIP4P/2005 model at T — and P = 0. The 
ice structures and variable labels are the same as those used in Tab. HI 



X (TIP4P/2005) 


Ih 


II 


III 


AX (Ih-II) AX 


(III-II) 


Vref (A 3 /molec.) 


31.34 


24.30 


25.26 


7.04 


0.96 


U s , re f (kJ/mol) 


-62.99 


-62.13 


-61.86 


-0.86 


0.27 


V (A 3 /molec.) 


33.20 


25.71 


26.76 


7.59 


1.05 


U s ,o (kJ/mol) 


-62.45 


-61.66 


-61.63 


-0.79 


0.03 


Uz,o (kJ/mol) 


16.15 


15.10 


16.29 


1.05 


1.19 


U (kJ/mol) 


-46.30 


-46.56 


-45.33 


0.26 


1.23 



Table III: Volume and energies of the studied ice phases as derived with the rigid TIP4PQ/2005 model at T — and P = 0. 
The ice structures and variable labels are the same as those used in Tab. [I] 



X (TIP4PQ/2005) 


Ih 


II 


III 


AX (Ih-II) AX (III-II) 


V re f (A 3 /molec.) 


30.67 


23.85 


24.92 


6.82 


1.07 


Us,ref (kJ/mol) 


-68.90 


-67.57 


-67.47 


-1.33 


0.10 


V (A 3 /molec.) 


32.51 


25.16 


26.32 


7.35 


1.16 


U s ,o (kJ/mol) 


-68.35 


-67.11 


-67.24 


-1.24 


-0.13 


Uz,o (kJ/mol) 


17.06 


15.91 


17.19 


1.15 


1.28 


U (kJ/mol) 


-51.30 


-51.20 


-50.05 


-0.10 


1.15 



tional contributions to the free energy are however sim- 
ilar. The coexistence pressure for the Ih-III transition 
at the temperature of T = 225 K is represented in Fig. 
[3] as a function of the relative static energy, AUs,ref, of 
the H-isomers of ice III. The relative energy is calculated 
with respect to the minimum potential energy of ice II, 
to allow for a comparison to available literature data. 7 
The coexistence pressure varies in the interval 0.18-0.24 
GPa at 225 K. There appears an approximate linear re- 
lation between the coexistence pressure and the static 
energy of the ice III isomer. The coexistence pressure 
reported using TI simulations deviates by less than 0.02 
GPa (about 8%) from the linear fit based upon the QH 
results. This small deviation suggests that the QHA is 
reasonably realistic even at this relatively high tempera- 
ture (T = 225 K). The difference between the QHA and 
the TI simulations is caused by the presence of anhar- 
monic effects not included in the QHA and also by the 
use of different H-isomers in both calculations. Unfortu- 
nately it is not possible to quantify the separate influence 
of both factors. 

An interesting result from our QH phase diagram is 
that, for cells with 324 molecules, ice II may be either 
stable at low T or unstable at all temperatures in the 
classical limit of the q-TIP4P /F water model. We find in 
Fig. [2] that ice II is unstable in the whole (T, P) region 
if the phase diagram is calculated with any of the three 
most stable H-isomers of ice III. Thus, the stability of 
ice II is determined by the static energy, Us.ref, of the 
H-isomer of ice III. 

Some differences in the phase diagrams calculated with 
TIP4P-like models might be caused by the differences in 



the static energy, Us, r ef-> of the single H-isomer chosen 
to represent ice III. For example, the classical phase di- 
agram for the rigid TIP4P model was calculated with 
a H-isomer of ice III with partial H-disorder^ We have 
seen in Fig. [T] that partial H-disorder is less stable (it 
has higher energy) than full H-disorder for TIP4P-likc 
models. Thus, this choice helps to increase the stability 
region of ice II. On the contrary, the phase diagram for 
the flexible q-TIP4P /F model was calculated with full H- 
disorder for ice III. Here the increased stabilization of ice 
III (see Fig. 1) plays an important role for the reported 
instability of ice II. 7 Our QH calculation strongly sug- 
gests that the ice II instability is not a deficiency of the 
q-TIP4P/F model, but a finite size effect related to the 
particular H-isomer randomly selected for the simulation. 

We turn now to the calculation of the QH phase di- 
agram of ice Ih, II, and III for the quantum limit of 
the q-TIP4P/F model. In this case there are quan- 
tum TI results for the melting of ice Ih at atmospheric 
pressure; 36 ' 37 but not for the coexistence between differ- 
ent ice phases. 



B. Quantum limit 

The QH phase diagram of ice Ih, II, and III in the 
quantum limit is plotted in Fig. 0J Finite size effects 
related to proton disorder are very important for ice III 
(N = 324). As in the classical limit, this effect is related 
to the differences in the static energy, Us,ref, of the H- 
isomers of ice III. Vibrational contributions to the free 
energy are however similar, and then the coexistence lines 
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Figure 3: Coexistence pressure for ice Ih and III at 225 K as 
a function of the relative static energy of the H-isomers of ice 
III. The static energy of ice II has been taken as zero of the 
energy scale. Circles are results derived by the QHA using the 
q-TIP4P/F model. The square is the result of Ref. Q based 
on the TI method with the same water potential. The line is 
a linear fit to the QH data. 



Figure 5: Comparison of the QH phase diagram of ice Ih, 
II, and III in the classical and quantum limits. Dotted lines 
correspond to experimental data from Ref. Q] that include the 
liquid (L) phase. Results derived with the q-TIP4P/F model. 
The calculation was performed with the H-isomer of ice III 
labeled with an arrow in Fig. [T] Circles show the position of 
the triple point for ice Ih-II-III. 
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Figure 4: QH phase diagram of ice Ih-II-III calculated with 
the q-TIP4P/F model in the quantum limit. The multiple 
coexistence lines II-III and Ih-III show the results for the six 
studied H-isomers of ice III. The phase boundary Ih-II is a plot 
of six superimposed lines, each one calculated with a different 
random H-isomer of ice Ih. Circles show the position of the 
triple point (TP) for the ice phases. 



calculated for the ice III isomers are nearly parallel. 

In contrast, the QHA reveals that finite size effects 
related to H-disorder are vanishingly small for ice Ih 
(N = 288). The coexistence line Ih-II was calculated 
for each of the six randomly generated H-isomers of ice 
Ih. In this case, the six calculated Ih-II phase boundaries 
appear superposed as a unique line at the scale of the fig- 
ure. The coexistence lines Ih-III have been displayed for 
a single H-isomer of ice Ih. 

Ice II is always a stable phase in the quantum phase di- 
agram independently of the employed H-isomer of ice III. 
A triple point Ih-II-III appears for all studied H-isomers, 
in contrast to the classical results in Fig. [5J The triple 
point temperature, Ttp, is found in an interval 75-136 
K depending upon the H-isomer of ice III. The triple 
point pressure, Ptp, appears in the interval 0.18-0.24 
GPa. The magnitude of these intervals provides a quan- 
titative estimation of the influence of the finite size effect 
of the proton disorder of ice III in the calculated phase 
diagram. We find that both quantities (Ttp, Ptp) show 
a linear dependence as a function of the static energy, 
Us, ref ^ °f i ce III, in a way very similar to that shown in 
Fig. [3] for the coexistence pressure between ice Ih and 
III. 
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Figure 6: Gibbs free energy difference of ice Ih and II as a 
function of P at T = K. Coexistence conditions are labeled 
by open symbols (AG = 0). Results derived with the q- 
TIP4P/F model in both classical and quantum limits. 



C. Comparison of quantum and classical limits 

The quantum and classical limits of the QH phase di- 
agram for the q-TIP4P/F model are compared in Fig. 
03 The H-isomer of ice III indicated by an arrow in Fig. 
[T] has been arbitrarily chosen for this comparison. The 
static energy of this H-isomer is the closest one to that 
of partial H-disorder structures. The most conspicuous 
quantum effect is the increased stability of ice II. The 
triple point Ih-II-III is found classically at (35 K, 0.3 
GPa), while the quantum limit is (136 K, 0.24 GPa), i.e. 
a shift of about 100 K and -0.06 GPa due to the consid- 
eration of quantum vibrational effects. The experimental 
boundaries in this region of the phase diagram are shown 
in Fig. [5] by dotted lines. The experimental triple point 
Ih-II-III is found at (239 K, 0.21 GPa)J- 



1. Coexistence Ih-II 

In the quantum limit, ice II occupies a large portion of 
the region of stability found classically for ice Ih. There- 
fore, the coexistence line Ih-II is shifted to lower pressures 
with respect to the classical one (see Fig. [5]). The Gibbs 
free energy difference between ice Ih and II, 

AG = Gih — Gij i (8) 

is plotted in Fig. [6] as a function of the pressure at T = 
K. The zero temperature condition implies that 

AG = AU + PAV . (9) 

The plot of AG in Fig. [S]is nearly linear in P. This fact, 
in the light of Eq. ©, implies that both AU and AV 
vary slowly with P in the studied pressure interval. 



Table IV: Average of the wavenumbers obtained with the q- 
TIP4P/F model for the volume Vo of the studied ices. Vo 
is the equilibrium volume in the quantum limit at T = 
and P = 0. The last two columns show the ratio of the 
wavenumbers with respect to the data of ice II. 
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3506 
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The main difference between the quantum and classical 
result for AG is the value of the ordinate at the origin, 
AUq. Fig. [5] shows that in the classical limit 



ADb.ela = - 1 - 14 kJ/mol (EE AU S ,ref) 



(10) 



a value that corresponds to the difference in the static 
energies, Us.ref, of ice Ih and II (see Tab. J]). In the 
quantum limit the ordinate at the origin is 



AU = -0.46 kJ/mol (= AU s ,o + MJ z ,o) 



(11) 



The quantum result differs from Allodia by an energy 
increment that essentially corresponds to the difference 
in the zero-point energy, Uz,o, of ice Ih and II at T = 
and P = 0. The data in Tab. [I] show that 



AU z ,o = Uz fi ,ih - Uzfiji = 0.67 kJ/mol 



(12) 



i.e., the zero-point energy of ice II is 0.67 kJ/mol lower 
than that of ice Ih. This is the physical reason for the 
quantum shift in the coexistence pressure of the Ih-II 
transition (abscissa of the open circles in Fig. [5]) and 
the origin of the increased stabilization of ice II in the 
quantum phase diagram. 

Why is the Uz,o of ice II lower than that of ice Ih? 
Uz,o is proportional to the average of the vibrational fre- 
quencies, aJfe, of the ice cell [see Eq. ©]. In Tab. IIV1 the 
average of translational, librational, bending and stretch- 
ing modes is presented for the equilibrium cells of the ice 
phases. We observe that the largest difference in oJk be- 
tween ice Ih and II is due to the librational modes, that 
are about 10% lower in ice II. The stretching modes show 
a competing behavior in the sense that they are larger for 
ice II. But the overall effect of all modes is the reduction 
of the zero-point energy of ice II in comparison to ice Ih 
by about 0.67 kJ/mol (1%). We can anticipate that for 
rigid water models this competing mechanism between 
librational and stretching modes is absent. Thus, the 
stabilization of ice II in the quantum phase diagram of 
rigid models should be even larger than in the flexible 
one. 



2. Coexistence II-III 

The Gibbs free energy difference between ice III and 
II is plotted as a function of the temperature in Fig. [7] 
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Figure 7: Gibbs free energy difference between ice III and II 
as a function of T at P = 0.33 GPa. Coexistence conditions 
are labeled by open symbols (AG = 0). Both classical and 
quantum limits were derived with the q-TIP4P/F model for 
the H-isomer of ice III labeled with an arrow in Fig. [1] 



Both classical and quantum QHA limits are displayed at 
a constant pressure of P = 0.33 GPa. The coexistence 
temperature in the quantum limit is shifted by about 
100 K toward higher temperatures, i.e., quantum effects 
play an important role in the stabilization of ice II. The 
zero-point energies in Tab. U show that Uz oi ice II is 
significantly lower (about 0.6 kJ/mol) than that of ice III. 
The slope of the AG curves in Fig. [7] is always negative 



d (AG) 
dT 



= -AS = Sti - Sttt <0 



(13) 



which is consistent with the larger entropy of ice III due 
to its H-disorder. At a given temperature the slope of 
the quantum AG curve is larger (in absolute value) than 
in the classical result. This implies that the excess of 
entropy of ice III, with respect to ice II, is larger in 
the quantum limit, i.e., the quantum vibrational entropy 
contributes to stabilize ice III. Nevertheless, the overall 
quantum effect in AG implies a strong stabilization of ice 
II with respect to ice III, as a consequence of its lower 
zero-point energy. 



V. RIGID TIP4P/2005 MODEL 

The QH phase diagram of ice Ih, II, and III with the 
rigid TIP4P/2005 potential has been calculated using the 
same H-isomers as those employed for the q-TIP4P/F 
study in Fig. Differences in the results should be 
attributed to the potential models (flexible versus rigid) 
and not to effects related to the selected H-isomers. 

The classical and quantum QH phase diagrams are pre- 
sented in Fig. [5] as dashed and continuous curves, respec- 
tively. The phase diagram reported for the TIP4P/2005 



Figure 8: Comparison of the QH phase diagram of ice Ih, 
II, and III in the classical and quantum limits. Dotted lines 
correspond to the classical TI simulations of Ref. that in- 
clude the liquid (L) phase. These results were derived with 
the TIP4P/2005 model. Ice III is modeled with the same H- 
isomer as that employed in Fig. [5] for the q-TIP4P/F model. 
Circles show the position of the triple point for ice Ih-II-III. 



model by classical TI simulations is shown by dotted 
lines. In the last case the coexistence with the liquid 
phase is also given. There are several aspects to be com- 
mented. First is the comparison between the classical 
QHA and the TI results. The coexistence lines between 
the ice phases are nearly parallel in both calculations. We 
recall that the slope of the coexistence lines is determined 
by the Clausius-Clapeyron relation 



dP coe AH 



dTr, 



AS 



T coe AV AV 



(14) 



where AH is the enthalpy difference (latent heat) be- 
tween the two phases at equilibrium. It has been shown 
earlier that the QHA provides a realistic approximation 
for the enthalpy and volume of ice Ih, II, and III in a 
broad (T, P) interval. 22 Therefore, the QHA shows a rea- 
sonable overall agreement to the classical TI results for 
the slopes of the coexistence lines. The largest devia- 
tion between the QH and TI slopes is found for the Ih-II 
transition at temperatures close to the triple point, where 
the slope of the QHA is lower. The larger stability re- 
gion of ice III in the QHA is likely due to the finite size 
effect related to the different H-ordering of the studied 
H-isomers, although the approximate treatment of an- 
harmonic effects by the QHA may be also the origin of 
this behavior. 

The consideration of quantum vibrational effects 
changes dramatically the QH phase diagram. The main 
quantum effect is the stabilization of ice II. It has two 
important consequences. The first is that ice II becomes 
the stable phase at low temperatures even at P = 0. 
The second is that ice III disappears as stable phase in 



the displayed region of the phase diagram. Both facts 
imply that the quantum QH phase diagram becomes in 
strong disagreement to experimental facts, as opposed to 
the classical one. 

In the case of the TIP4P/2005 model there are not 
quantum TI results available for comparison. Neverthe- 
less the QH prediction that ice II becomes the stable low 
temperature phase agrees with the quantum path integral 
simulations of ice Ih and II in Ref. [||. The explanation 
given for this behavior was that the TIP4P/2005 model 
is parameterized to be used in a classical limit. The com- 
bination with quantum simulations implies some kind of 
overcounting of quantum effects that leads to unphysical 
results. 

It is interesting to analyze the physical reason for the 
larger stabilization of ice II in the quantum phase di- 
agram of Fig. [8] (rigid model) in comparison to Fig. 
|S] (flexible model). At T = and P = the shift 
in the coexistence pressure of ice Ih-II in the classical 
and quantum limits is proportional to the difference in 
the zero-point energy, AUz,o, between both phases. For 
TIP4P/2005 we get (see Tab. ED) 

AU z .o = Uz fi jh - Uz,o,n = 1-05 kJ/mol , (15) 

i.e., the zero-point energy of ice II is more than 1 kJ/mol 
lower than that of ice Ih. This stabilization is about 
50% larger than that of the flexible model [see Eq. ([12"]) ]. 
This fact is a consequence of the absence in a rigid model 
of competing contributions of librational and stretching 
modes to AUzo, as discussed in Subsec. IIV C II In con- 
clusion, the stabilization of ice II by its lower zero-point 
energy is larger for the rigid model than for the flexible 
one. As a consequence ice II becomes the stable phase at 
low temperature and ice III is unstable in the quantum 
limit of the TIP4P/2005 model (see Fig. 0). 

A last remark on the comparison of the QH phase dia- 
grams of the rigid and flexible models. The most realistic 
results are obtained in the quantum limit for the flexible 
model (see Fig. [5j , but in the classical limit for the rigid 
model (see Fig. [5]). The models were parameterized to 
be used in either quantum (the flexible one)^ or classical 
simulations (the rigid one)^ Thus, the best result corre- 
lates in each case with the conditions where the model 
was parameterized. 



VI. RIGID TIP4PQ/2005 MODEL 

The rigid TIP4PQ/2005 model differs from 
TIP4P/2005 by an increase of about 4% in the 
point charges, i.e., the parameter qn is the only one 
that changes from 0.5564 e to 0.5764 e* The QH phase 
diagram of ice Ih, II, and III for the rigid TIP4PQ/2005 
potential has been calculated with the same H-isomers as 
those used in the studies shown in Fig. [SJ (q-TIP4P/F) 
and Fig. ]5](TIP4P/2005). The results for TIP4PQ/2005 
are displayed in Fig. [5] 



Oh 



— i 1 1 1 r-n 1 1 1 1 1 1 — 

'^il-lll 1 
■ " Ih-II Ih-III^---- 


- 
- 






- 


classical QHA 

quantum QHA 

o TP (Ih-II-III) 




1 1 1 1 1 1 1 


TIP4PQ/2005 
j i 1 i i i 





100 200 300 

T(K) 



Figure 9: Comparison of the QHA phase diagram of ice Ih, II, 
and III in the classical and quantum limits. Results derived 
with the TIP4PQ/2005 model. Ice III is modeled with the 
same H-isomer as that employed in Fig. [8]for the TIP4P/2005 
model. A circle shows the position of the triple point for ice 
Ih-II-III in the classical limit. 

We observe that the QH diagrams in Figs. [8] and [9] 
are qualitatively similar. In the classical case, there ap- 
pears a triple point for the ices Ih-II-III. However, in the 
quantum limit, ice III is unstable and the phase diagram 
displays only the ice Ih-II transition. Another similarity 
of both phase diagrams is that the coexistence lines in 
Figs. [5] and [5] are nearly parallel. 

The change in the point charges modifies the static en- 
ergy, Us.ref, of the ice phases (see Tabs. ITlland Ull]) . This 
translates into a shift of the coexistence lines calculated 
with both rigid models. For example, in the quantum 
limit the transition Ih-II appears at higher pressures for 
TIP4PQ/2005, so that ice Ih becomes the stable low tem- 
perature phase (see Fig. |9]). This result is in agreement 
with the path integral simulations of the TIP4PQ/2005 
model in Ref. [^. 

Classical and quantum phase diagrams of water, de- 
rived by TI simulations with the TIP4PQ/2005 model, 
have been reported recently^ The results differ markedly 
from our QHA. We believe that the main reason of dis- 
crepancy is the finite size effect in the value of the static 
energy, Us. ref, of ice III. As a check of this hypothesis, we 
have repeated our QHA calculation of the TIP4PQ/2005 
phase diagram by shifting the internal energy of ice III by 
a constant amount. All other model parameters remain 
unchanged (i.e., static energies of ice Ih and II, and the 
vibrational properties of the three ice phases). We have 
analyzed the effect of setting the static energy of ice III 
as 

U S ,ne W {V) = U S (V) + AU ref , (16) 




Figure 10: Comparison of the QH phase diagram of ice Ih, II, 
and III in the classical and quantum limits. Results derived 
with the TIP4PQ/2005 model. The results differ from those 
in Fig. [9] by an additional stabilization of ice III by a constant 
energy shift given by AU re f- 

where AU re f is a constant energy shift. The QH phase 
diagram for AU re f = was displayed in Fig. [9] The 
phase diagrams derived for AU re f — —0.29 kJ/mol and 
AU re f = —0.47 kJ/mol are displayed in Figs. [TU1 and 
ITTT respectively. The differences between these phase di- 
agrams are caused by the artificial stabilization of ice III 
by the constant energy AU re f. Note that the selected 
energy shifts are in the order of the dispersion of Us,ref 
represented in the abscissa of Fig. [Tjfor the q-TIP4P/F 
model. 

In the quantum limit, the phase diagram calculated 
with AU re f = -0.29 kJ/mol (Fig. [JO), is identical to 
that shown for AU re f — (Fig. [9]). The coexistence Ih-II 
is the only transition. The additional stabilization of ice 
III is not large enough to make it stable in the quantum 
limit. However, it does change drastically the classical 
limit of the phase diagram. Ice III becomes more stable 
than ice II in the whole region, and the classical phase di- 
agram loses its triple point and now displays only the ice 
Ih-III transition. Note that for AU re f = -0.29 kJ/mol, 
the triple point Ih-II-III is missing in both classical and 
quantum limits. 

Imposing a larger stabilization to ice III (AU re f — 
—0.47 kJ/mol), we obtain the QHA phase diagram shown 
in Fig. QT] In this plot we have represented also the clas- 
sical and quantum phase diagrams for ice Ih, II, and III 
calculated by TI simulations in Ref. |g. We observe that 
this additional stabilization brings the QHA in reason- 
able agreement to the TI results reported for this model. 
Now a triple point Ih-II-III is observed in the quantum 
limit, while the transition Ih-III is the only line in the 
classical case. Furthermore, we note that the slopes of 



Figure 11: Comparison of the QH phase diagram of ice Ih, II, 
and III in the classical and quantum limits. Results derived 
with the TIP4PQ/2005 model. The results differ from those 
in Fig. [9] by an additional stabilization of ice III by a constant 
energy shift given by AU re f- Dotted (dashed-dotted) lines 
correspond to the classical (quantum) TI simulations of Ref. 
6. Circles show the position of the triple point for ice Ih-II-III. 

the coexistence lines predicted by the QHA are in rea- 
sonable agreement to the TI results. 

The sequence of phase diagrams shown in Figs. |9j 
I10[ and [TJJ provides a vivid illustration of the dramatic 
changes in the ice phase diagram as a function of the 
stability of the employed H-isomer of ice III. 

VII. AVERAGE OVER PROTON DISORDER 

In this Section we comment on the necessity of per- 
forming some form of proton disorder average, at least 
for ice III, prior to the calculation of the phase diagram. 
Assuming that the number of H-isomers for an ice cell is 
M, the canonical partition function of the ice phase can 
be expressed as 

M 

e -?F = Y j e-^ Us -+ F ^ , (17) 

i=l 

where Us,i and F v i are the static energy and the vibra- 
tional free energy of i'th H-isomer. Note that many of 
the M H-isomers might be energetically degenerate as 
a consequence of the lattice symmetry. Eq. (|17[) is the 
formally correct way to average over the proton disorder 
of the ice phase. It has been applied to study order- 
disorder transitions of ice phases using small units cells, 3 
Obviously for large unit cells the total number M of H- 
isomers grows in such a way that the application of Eq. 
(fTT)) becomes an impossible task. 
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The alternative for large unit cells is the use of Eq. ([IJ , 
which is applied to a single H-isomer selected randomly 
from the set of M available ones4^— Note that in this 
equation the average over proton disorder is introduced 
ad hoc by the term with the proton disorder entropy Sh- 
An implicit assumption in Eq. (JlJ is that the cell is so 
large that the static energy, Us, and the vibrational free 
energy, F v , of the single H-isomer do not require any 
further average over the proton disorder. This assump- 
tion is correct in the thermodynamic limit as the relative 
fluctuation of thermodynamic quantities is expected to 
decrease as 1/y/N. 

However, we have shown that, for typical cell sizes used 
in simulations, the fluctuation of Us for ice III is far from 
its ideal thermodynamic limit. Therefore some form of 
proton disorder averaging of Us is necessary. A com- 
putational feasible proposal is to average only the term 
that shows the largest fluctuation as a function of the 
H-disorder, which is the potential energy of the refer- 
ence cell Us, re f associated to the chosen H-isomer. Thus, 
a simple proposal to average over proton disorder is to 
modify the static energy in Eq. ([I} by 

U S ,ave{V) = U S (V) + AU ave . (18) 

AU ave here is a constant energy shift that modifies the 
stability of the single selected H-isomer of ice III by an 
amount determined by the average U s,ref calculated over 
a random set of H- isomers, i.e., 

AUave = Us,ref — Us,ref ■ (19) 

Note that the proposed disorder averaging can be applied 
to the calculation of phase diagrams either by the QHA 
or by TI simulations. 

The average of Us,ref over the set of six random H- 
isomcrs of ice III with full H-disorder studied in this work 
gives Us.ref = -60.96 ± 0.04 kJ/mol for the q-TIP4P/F 
model. We estimate that the error of the mean value, 
Us,ref, should be as low as 0.01 kJ/mol for a reasonable 
convergence over the proton disorder. Then the size of 
our random sampling should be increased by a factor of 
16 to reduce our estimated error to this limit, i.e., for 
an ice III cell with 324 molecules one should increase 
the sampling of Us,ref to about 100 H-isomers. By using 
larger units cells, e.g., a 4 x 4 x 4 supercell with 768 water 
molecules, the number of required H-isomers to obtained 
a converged value of U s,ref should be lower than this. 
However, in terms of computational efficiency, the lower 
number of H-isomers may be overcompensated by the 
higher computational cost in the minimization of the cell 
energies. 

In the case of the partially disordered ice III we get 
U s,ref = — 60.73 ± 0.01 kJ/mol for our set of six random 
H-isomers with the q-TIP4P/F model. This value has 
achieved already the desired convergence. A last com- 
ment is that the proton disorder entropy Sh for par- 
tially disordered phases is lower than the Pauling esti- 
mate in Eq. @. The estimation of Sh for the frac- 
tional H-occupancies experimentally determined for ice 



III amounts to about 90% of the Pauling result. This 
entropy lowering has a significant influence in the phase 
diagram^ Numerical methods for the estimation of the 
proton disorder entropy of partially H-disordered phases 
can be found in Refs. [35U381 



VIII. CONCLUSIONS 

In this work we have presented a detailed study of the 
phase diagram of ice Ih, II, and III calculated with the 
QHA and TIP4P-like models. Several advantages of the 
QHA are worth to be mentioned: its computational cost 
is low, it is free from the statistical errors inherent to 
any numerical simulation, it can be applied to both clas- 
sical and quantum limits, and the most accurate results 
are expected in the low temperature limit, where anhar- 
monicities are lower. These advantages make the QHA an 
appropriate option to study finite size effects that can be- 
come prohibitively expensive in numerical (Monte Carlo 
or molecular dynamics) simulations. 

The effect of proton disorder in the phase diagram of 
TIP4P-like models has been a focus of our study. We 
have found that for the typical cell sizes employed in 
computer simulations, the finite size effect of H-disorder 
in ice Ih is very small. However, this effect is very large 
for ice III, so that the transitions II-III and Ih-III are 
strongly affected by it. The physical reason for this be- 
havior is that the static energy of ice III may change by 
an amount of several tenths of k J /mol depending on the 
considered H-configuration. Thus a randomly selected ice 
III structure makes the calculated phase diagram affected 
by an uncontrolled factor that can be highly significant 
for the final result. Crystallographic data of ice Ih and 
III reveal significant differences as a consequence of the 
larger structural complexity of ice III. As example, ice Ih 
is a network of hexagon rings of O-atoms, while in ice III 
there appear five-, seven, and eight-members rings. The 
tetragonal 0-0-0 angles deform from the ideal value of 
about 109° in ice Ih into angles between 80° and 140° 
in ice III.— Thus the large finite size effect related to 
the H-disorder in ice III, in comparison to ice Ih, seems 
to correlate with its increased structural complexity. We 
have discussed the necessity of performing a disorder av- 
erage of the lattice energy in order to reduce this effect. 

Another aspect related to the H-disorder in ice III is 
that TIP4P-likc models predict that full H-disorder is 
more stable than partial H-disorder. However, this result 
is against the data derived from diffraction experiments 
of ice Ill^i A phase diagram derived with ice III hav- 
ing partial H-disorder may be significantly different from 
that derived with full H-disorder, as a consequence of the 
differences in their static energies. 

These findings allows us to rationalize previously con- 
tradictory results of phase diagrams calculated with 
TIP4P-like potentials. A significant example was the 
reported instability of ice II in the classical limit of 
the flexible q-TIP4P/F model. 7 This fact contrasts 
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with the stability of ice II reported with rigid TIP4P 
and TIP4P/2005 models^ Our QH results show that 
the stability of ice II strongly depends upon the H- 
configuration of the chosen ice III isomer. The fact that 
ice III has been described either with partial 4 or full H- 
disorder— i makes difficult the comparison of phase dia- 
grams with different models. It is not easy to discrimi- 
nate from the reported simulations the effects due to the 
variable stability of the ice III isomers from those caused 
by the different water models. 

The QH phase diagram of ices Ih-II-III for the q- 
TIP4P/F, TIP4P/2005 and TIP4PQ/2005 effective po- 
tentials has been able to reproduce qualitatively most 
features that had been previously studied by TI in both 
classical and quantum simulations. Our results are ob- 
tained using the same H-isomers for the three TIP4P-likc 
models, allowing for an easier interpretation of the differ- 
ences encountered in the phase diagrams. We have found 
that for the flexible model (q-TIP4P/F) the triple point 
Ih-II-III is shifted in the quantum limit by 100 K and 
-0.06 GPa with respect to the classical one. This effect 
is very large, specially in the temperature. Its physical 
origin is related to the lower zero-point energy of ice II, 
when compared to that of ice Ih and III. This fact trans- 
lates in an increased stability of ice II when vibrational 
quantum effects are considered. The average frequency 
of molecular librations in the H-bond network are nearly 
10% smaller in ice II than in the other ice phases. This 
causes a significant reduction of the zero-point energy. 
Interestingly, the intramolecular stretching modes of ice 
II are predicted at higher frequencies than those in ice 
Ih and III. This anticorrelation between libration and 
stretching modes has been often stressed in the litera- 
ture, i.e., any factor that shifts librational frequencies 
in one direction acts also modifying the stretching fre- 
quencies in the opposite one^ We find as net effect that 
librational modes dominate over stretching ones, leading 
to the stabilization of ice II by its lower zero-point energy. 

The main difference between the QH phase diagrams 



obtained by the flexible and rigid models is associated 
to the absence of anticorrelation between H-bond libra- 
tions and O-H stretchings in the rigid models. Obvi- 
ously intramolecular bonds are frozen for rigid water. 
As a consequence, the stabilization of ice II by its zero- 
point energy is larger for the rigid models (TIP4P/2005, 
TIP4PQ/2005) than for the flexible one (q-TIP4P/F). It 
may be even so large that ice II becomes the stable phase 
at low temperatures. This result has been reported by 
quantum path integral simulations with the TIP4P/2005 
potential in Ref. 5 and is also reproduced by our QHA. 
Another effect related to the large stabilization of ice II is 
that this phase may occupy the stability region of ice III 
in the quantum limit of the rigid models. Then the triple 
point Ih-II-III does not appear in the quantum phase 
diagram. This unphysical behavior, as displayed in the 
quantum phase diagrams of Figs. [8]and[9l can be avoided 
if the H-isomer of ice III is particularly stable, as was 
shown in the phase diagram of Fig. [TT1 

The present work can be extended along several di- 
rections. An obvious one is the analysis of the phase 
diagram of other ice phases with the flexible q-TIP4P/F 
model as well as the study of finite size effects in the 
proton disorder of other phases, as ice V, VI, and VII. 
A second aim is the use of DFT to avoid the limitations 
of the empirical potentials, in particular with respect to 
the energetics associated to proton rearrangements in ice 
phases. 
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